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摘 要 椭圆 函数 是 一 种 特殊 的 双 周 期 复 变 函数 ， 广 泛 应 用 于 工程 问题 中 ， 尤 其 非 线性 问题 中 居 
多 .在 工程 中 遇 到 的 椭圆 函数 以 二 阶 椭圆 函数 为 主 ， 而 且 很 多 复杂 的 椭圆 函数 都 可 以 通过 变换 由 二 
阶 椭圆 函数 得 到 .二 阶 椭圆 函数 包括 Jacobi 椭圆 函数 和 Weierstrass 椭圆 函数 .它们 都 可 以 进行 罕 
级 数 展开 ， 直 接 计 算 很 不 方便 ， 椭 圆 函 数 的 重要 性 质 之 一 就 是 具有 加 法 定理 ， 所 以 可 以 利用 精细 积 
分 法 0 求解 .通过 使 用 精细 积分 法 成 功 求解 了 二 阶 椭 圆 函 数 ， 并 指出 精细 积分 法 是 一 种 求解 椭圆 
函数 的 高 效 稳定 的 方法 . 


关键 词 ”精细 积分 ， Jacobi 椭圆 吕 数 ， Weierstrass 撕 贺 马 数 ， 双 周期 
引言 


精细 积分 法 首先 解决 了 时 不 变 系 统 的 时 间 积 分 问题 ， 并 取得 了 达到 计算 机 精度 的 高 精度 结 
果 . 精细 积 分 法 有 2 个 要 点 (1) 运用 2 类 算法 有 ?9; (2) 把 注意 力 放 在 增 量 上 ， 而 不 是 全 
量 ， 解 决 时 不 变 系统 之 后 ， 各 种 近似 方法 被 广泛 的 应 用 于 求解 其 他 问题 ， 例 如 时 变 系 统 或 非 线 
性 系统 的 时 间 积分 等 . 

在 工程 问题 中 常 遇 到 椭 俩 函数 数值 求解 的 问题 ， 该 问题 可 以 通过 查 表 和 使 用 软件 来 解决 ， 
例如 Matlab 和 IMSL 等 软件 就 有 求解 椭圆 函数 的 模块 ， 本 文 利用 精细 积分 法 实现 了 对 Jacobi 
椭圆 函数 和 Weierstrass 椭圆 函数 的 数值 求解 ,大 量 的 数值 算 例 指出 精细 积分 法 得 出 的 结果 具有 
极 高 的 精度 ， 而 且 无 论 是 在 效率 上 还 是 在 稳定 性 上 精细 积分 法 都 明显 的 优 于 IMSL 数学 库 的 椭 
圆 函数 模块 . 


1 Jacobi 椭圆 函数 的 精细 积分 算法 
第 二 类 Legendre 椭圆 积分 可 以 定义 为 


3 dz 
2 | (1 — z2)(1 — kz2) . 


Jacobi 椭圆 函数 snu = z 定义 为 第 二 类 Legendre 椭圆 积分 的 反 苍 数 ， 男 外 两 个 基本 Jacobi 椭 
圆 函数 分 别 定义 为 


cnu= Vl-— sn? = V1— 22, dnu = V1 — k?sn2u = TR (2) 
其 中 ， 参 数 k 称 为 模 数 ， 这 三 个 Jacobi 椭圆 函数 的 寡 级 数 展 开 式 分 别 为 外 
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1 
snu = (+) + (+ 14k? + kt)us — 
1 2 1 2 4 
nu=1- + 广 (1+ 和 多 Ja — Bi(l + 44k + 16k4)us 十 .… (3) 


dnu=1— kt 十 一 本 (4 十 1 - (6k? + 44k4 十 ke)ws 十 .…- 


相应 的 加 法 公式 为 
sn(2v) = 2snv cn : dnv cn(2v) = cm au 一 sm dmn A dmn2u — kn2v . sn2v .cn2v 


1 一 K2 .8s724V 1—k2.sntiv 1 — k2. sn4v 


(4) 

现在 需要 解决 的 问题 是 ， 给 出 坐标 变量 w 和 模 数 k 数值 求解 这 三 个 Jacobi 椭圆 函数 ， 利 用 这 
些 函 数 的 加 法 定理 (4), 精细 积分 法 可 以 用 于 求解 该 类 问题 . 

为 了 进行 数值 积分 ， 引 入 一 个 小 的 空间 步 长 ， 记 为 : 7. 于 是 产生 一 系列 等 步 长 的 空间 间 


隅 为 
wo = 0， ul 三))， Ug = Kn, "1 U2L =2rn=u (5) 


假设 这 三 个 Jacobi 椭圆 函数 位 于 ”点 的 值 都 已 经 求 得 ， 那 么 通过 利用 加 法 定理 (4) 就 可 方便 的 
求 出 函数 位 于 点 的 值 ， 因 此 ， 原 问题 变 为 在 ” 点 对 Jacobi 椭圆 函数 进行 数值 求解 ， 令 


T= 二 n/m, 其 中 m= 2w, 例如 N=20，m = 1048576 (6) 


由 于 7 本 来 是 不 大 的 空间 区 段 ， 则 7 = n/m 将 是 非常 小 的 一 个 空间 区 段 了 . 因此 对 于 空间 区 段 
7T, 有 


SNT AST 一 zi 十 上 2)73 十 (1 +14k2 二 k4)7? 一 9; 


1 1 
cnr ws 一 可 72 不 (1 + 4k?)74 -+44k +16k)r =1+C! (7) 


dnT 守 1— 于 


(4 + k4)T4 一 (16k? + 44k4 二 16)rs =1+Di 

因 7 很 小 ， 帮 级 数 4 项 展开 应 以 足够 .在 上 式 中 ， Su C1 和 D1 相对 于 单位 1 是 非常 小 的 量 . 

因此 在 计算 过 程 中 至 关 重 要 的 一 点 是 数值 的 存储 只 能 是 (7) 式 中 的 51, C1 和 Di, 而 不 是 1+ Ca 

和 1+ Di. 因为 C1 和 Di 很 小 ， 当 它们 与 单位 1 相 加 时 ， 就 会 成 为 尾数 ， 在 计算 机 的 舍 入 操作 

中 ， 其 精度 将 丧失 殖 尽 . S51, C1 和 Di 就 是 增 量 “入 . 这 就 是 以 上 提 到 的 精细 积分 法 的 第 二 个 
把 (7) 式 代 入 (4) 式 ， 并 令 snv = 5; 可 以 得 到 


2 


= 92; 


k2*St +2C;+ OC? — SS?(1+ Di) 


cn(2v) =1+ Fc{(Ci)=1+ 1 — k253 


= 1+ C2; (8) 


KS +2Di+ D? — k2S?(1 + Ci) 


dn(2v) = 1 + Fp(D:;) = 1 R25 


= 1+ D2: 
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因此 ， 通 过 (8) 式 我 们 可 以 得 到 Sm, Cm 和 Dm, 从 而 求 得 sn(n), cn(n) 和 dn(n). 这样 的 计算 一 
共 需 要 执行 N 次 ， 而 且 只 有 Si, Ci 和 Pi (i = 2 ,2 ,2”) 被 计算 并 存 入 内 存 . (8) 式 的 入 
次 计算 相当 于 一 下 语句 


for (iter = 0; iter < Ni iter + + )i=2iter, 52; = Fs(S$i); C2; = Fc(Ci;), Do = Fp(D;) 
(9) 
当 以 上 语句 的 循环 结束 后 ， 再 执行 


sn(m7T) = sn(n) = Sm, cn(m7T)= en(n) =1+Cn, dn(m7) = dn(n)=1+ Dn (10) 


便 可 .由 于 N 次 登 代 后 ， Sm, Cm 和 Dm 已 不 再 是 很 小 的 量 了 ， 这 个 加 法 已 没有 严重 的 合 入 误 
差 了 电光 . 以 上 便 是 Jacobi 椭圆 函数 的 精细 积分 算法 . 


2 Weierstrass 椭圆 函数 的 精细 积分 算法 
第 一 类 Weierstrass 椭圆 积分 定义 为 


=/ °° (11) 
oo V4z — g2z— g3 


Weierstrass 椭圆 函数 z = p(w) 定义 为 第 一 类 Weierstrass 椭圆 积 分 的 反 项 数 . 上 式 中 的 gz, gs 
被 称 为 Weierstrass 椭圆 函数 的 不 变量 ,它们 是 Weierstrass 椭圆 函数 的 2 个 周期 ww 和 十 的 函数 


， 140 ;140 
92=6054= 2 Ti, 9=14056= > (12) 


mm! m,m! 


其 中 ， w= 2mw 十 2m/w'; 式 (12) 中 的 mm 和 m' 是 对 全 体 整 数 求 和 ,但 m 和 m 不 能 同时 为 0. 
Weierstrass 椭圆 函数 的 精细 积分 算法 与 第 一 节 中 得 到 的 Jacobi 椭圆 防 数 的 精细 积分 算法 
相 类 似 . 因此， 这 里 仅 阐 述 该 算法 的 主要 过 程 与 公式 . 
Weierstrass 椭圆 函数 的 寡 级 数 展开 式 和 加 法 定理 分 别 为 巴 


1 9 9 , 9 6  g293 8 
i eA 23 22253 ye 13 
gt) 一 五 十 20 二 28 十 12002 + 560" 十 (9 


PA(vV) 十 39g2p2(v) 十 2gsp(u) 十 让 9 


本 493(v) — gap(v) — gs s 
给 出 一 个 极 小 的 空间 区 段 ， 7 = 7/m, 其 中 m= 2 ,并 令 
2 
P= 7? 十 i 十 1 十 er (15) 
于 是 ， 方 程 (13) 可 以 重 写 为 ; 
p(T) OT -s+h (16) 


把 (16) 式 代入 (14) 式 ， 并 令 ， p(v) = 元 上 + 有 ,= 2i7, 加 法 定理 (14) 可 以 重 写 为 


3 9 1 
1 1 -ga(v + vSP)? + -gv + vi Pi) + mg2v" 
FD, RD 
(2v)? 4 a1 + viP) — gv + voP) — vegs 
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+ F(v,P:) = + Poi (17) 


1 
(2v)? (2o)? 
在 式 (17) 中 已; 是 常规 部 分 。 Weierstrass 椭圆 函数 的 精细 积分 算法 可 以 校 下 面 过 程 实现 
首先 ， 按 照 以 下 语句 进行 N 次 循环 计算 


for (iter = 0; iter < N; iter + 二 +)i=2", Pi= F(i7,P) (18) 


注意 只 有 常规 项 已 ; 的 值 而 不 是 pli7) 的 值 存 入 计算 机 内 存 . 
然后 ， 当 进行 完 以 上 的 循环 计算 时 ， Weierstrass 李峰 蚂 数 的 数值 结果 可 以 由 下 式 得 到 


p(w) = p(mr) = Pn + 广 (19) 


因此 ， 通 过 以 上 推导 我 们 得 到 了 Weierstrass 椭圆 函数 的 精细 积分 算法 . 


3 数值 算 例 与 误差 分 析 


在 本 节 中 将 讨论 精细 积分 算法 的 精度 , 效率 和 稳定 性 ， 相应 的 误差 分 析 也 将 给 出 . 在 本 节 中 
精细 积分 算法 的 程序 采用 Fortran90 语言 编写 , 而且 在 全 部 的 算 例 中 如 无 特殊 说 明 都 取 N = 20. 


3.1 Jacobi 椭 加 函数 
(A) 晓 化 为 三 角 函 数 演 算 
当 Jacobi 椭圆 函数 的 模 数 & = 0 时 ， Jacobi 椭圆 函数 就 晓 化 为 三 角 函 数 


snulkg=0 = sinu, cnulk=0= cosu, dnulk=0 = tanu (20) 


三 角 鲨 数 的 精确 数值 求解 早已 在 各 软件 和 计算 机 语言 中 实现 ， 通 过 和 相应 的 三 角 函 数值 对 比 ， 
就 可 以 得 出 精细 积分 法 相对 精度 . 

在 复 域内 取 500000 个 点 (不 包括 极点 ), 用 精细 积分 法 计算 Jacobi 椭圆 函数 当 = 0 时 在 
各 点 的 值 ， 并 用 各 点 相应 的 三 角 函 数值 对 比 ， 结 果 给 出 全 部 点 的 相对 误差 均 小 于 10-. 

(B) 特殊 点 验算 

Jacobi 椭圆 星 数 在 每 个 单 胞 内 都 存在 特殊 点 


snK =1, oe dnK =k’ 

f 
sn(K +iK')= 下 ， cn(K +ikK')= 1 dn(K +iK')= (21) 
sn(iK’') = oo， i 一 co， dn(iK') = 


上 式 中 的 kW 称 作 补 模 数 ， K 和 K' 为 分 别 与 上 和 ww 相对 应 的 全 椭圆 积分 ， 相 关 的 定义 可 以 参 
看 文 四 

通过 对 取 不 同 模 数 的 500000 点 的 验算 ， 发 现 相对 误差 的 变化 范围 为 10-6 ~ 10-16. 分 
析 相对 误差 产生 的 原因 发 现 是 由 于 ， 当 大 和 k' 接近 0 或 1 时 ，K 和 K' 无 法 得 到 精度 高 的 近 
似 解 

(C) 恒等式 验算 

Jacobi 椭 立 了 艺 数 存 在 如 下 恒等式 


sniut+cniu=1, dniut+kisn2u=1 (22) 
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对 于 不 同 模 数 的 Jacobi 椭圆 函数 在 复 空间 内 的 任意 一 点 (不 包括 极点 ) 都 应 该 满足 (22) 式 . 
在 一 个 单 胞 内 取 500 000 个 不 同 的 坐标 点 作为 算 例 ， 计 算 Jacobi 椭圆 函数 的 值 并 回 代 (22) 式 得 
到 ， 所 有 算 例 的 相对 误差 均 小 于 10-14. 
3.2 Weierstrass 椭圆 函数 

Weierstrass 椭圆 函数 没有 类 仅 于 Jacobi 椭圆 梢 数 的 恒等式 ;而 且 对 于 给 定 的 9a, 9s 的 情 
况 无 法 求 得 该 函数 的 周期 2w 和 2w' 所 以 直接 的 数值 检验 方法 较 少 .Weierstrass 椭圆 函数 与 
Jacobi 椭圆 项 数 之 间 存 在 如 下 关系 


sn2v = 人 Re (23) 
v el 一 e2 
o( -十 = =) €2 
式 中 的 ei,e2 和 es 称 作 Weierstrass 椭圆 函数 的 平稳 值 ， 同 时 也 是 下 面 三 次 方程 得 根 
4z — gaz—g3=0 (24) 


由 于 Jacobi 椭圆 函数 已 经 可 以 精确 求解 ， 故 可 以 利用 式 (23), 式 (24) 式 将 求 得 的 Weierstrass 椭 
圆 函 数 与 相应 的 Jacobi 椭圆 孙 数 对 比 ,进而 验证 Weierstrass 椭圆 函数 的 精度 . 同样 计算 500 000 
个 算 例 ， 得 到 ， 所 有 算 例 的 相对 误差 均 小 于 10-14. 

对 于 给 定 周期 2w 和 2w' 的 情况 (采用 的 精细 积 算法 与 第 2 节 中 所 论述 的 不 同 ， 由 于 篇 幅 
关系 就 不 在 这 里 详细 论述 ， 具 体 算法 将 在 另 一 篇 文章 中 给 出 ), 就 可 以 进行 周期 性 验证 ， 在 复 空 
间 内 取 一 对 坐标 (71, 胃 1) 和 (x2,92), 这 两 个 坐标 点 的 空间 间隔 为 


(z2 + yat) — (zi1 + 7) = 2 + 2:2 (25) 


上 式 中 4 和 ls 为 整数 ， 且 全 | + 人 | = 10. 根据 椭圆 隆 数 的 双 周 期 性 ， 孙 数 应 该 在 这 两 点 取 相 
同 的 项 数值 ， 分 别 计 算 两 点 的 值 并 对 比 就 可 以 得 到 相应 的 相对 误差 . 
验算 500000 对 坐标 点 ， 得 到 相对 误差 均 小 于 10-…. 


3.3 与 IMSL 计算 模块 对 比 

现在 很 多 软件 都 存在 计算 Jacobi 椭圆 函数 的 模块 ， 例 如 Matlab,Mathematica 和 IMSL. 这 
里 我 们 将 Jacobi 椭圆 消 数 的 精细 积分 法 程序 与 IMSL fortran 库 的 Jacobi 函数 模块 在 进行 一 下 
对 比 . 
3.3.1 计算 速度 对 比 

分 别 用 IMSL 和 精细 积分 法 将 200 000 个 算 例 计算 3 次 ， 并 将 每 种 方法 所 消耗 的 时 间 计 入 
表 1. 


表 1 IMSL 和 精细 积分 的 效率 对 比 


消耗 的 时 间 (s) 
IMSL 精细 积分 法 
N=10 N = 20 
第 1 次 4.13600000000000 1.68300000000000 3.05400000000000 
第 2 次 4.05500000000000 1.75200000000000 3.09400000000000 
第 3 次 4.11600000000000 1.76300000000000 3.05500000000000 


当 N = 20 时 ， IMSL 和 精细 积分 法 的 精度 都 可 以 达到 双 精 度 的 要 求 ， 但 从 表 1 中 可 以 看 
出 ， 精 细 积 分 法 的 效率 大 约 是 IMSL 的 1.3 倍 ， 而 且 调 用 精细 积分 法 求解 时 ， 一 次 计算 可 以 同 
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时 给 出 sn,cn,dn 的 值 ， 这 就 意味 着 如 果 一 个 包含 椭圆 函数 的 表达 式 中 同时 需要 求解 同一 点 的 
不 同 Jacobi 椭圆 函数 值 那 么 精细 积分 法 的 效率 就 可 以 在 提高 2 倍 (同时 计算 2 种 Jacobi 椭圆 
函数 ) 或 3 倍 (同时 计算 3 种 Jacobi 椭圆 函数 ), 此 时 精细 积分 法 的 效率 就 是 IMSL 的 2.6 或 3.9 
信 了 ， 这 是 非常 可 观 的 ， 如 果 该 问题 需要 大 量 达 代 则 可 以 节约 很 多 时 间 . 
3.3.2 稳定 性 对 比 

Jacobi 椭圆 函数 在 每 个 单 胞 内 存在 2 个 奇 点 ， 在 靠近 奇 点 的 区 域 是 无 论 是 IMSL 模块 还 是 
精细 积分 法 都 会 产生 误差 . 如 果 在 复 空间 划分 一 个 区 域 : (10 一 10i,10 一 10i, 10 十 10i, 一 10 十 102)， 
相对 模 数 有 = 0.5, 在 该 区 域内 3 种 Jacobi 椭圆 函数 都 存在 完整 的 单 胞 .在 该 区 域内 沿 实 轴 等 
间距 分 别 作 n -1 条 虚 轴 的 平行 线 ， 同 理 沿 虚 轴 也 作 n 一 1 条 实 轴 的 平行 线 ， 空间 被 刨 分 出 
(+1l) x+1l 个 点 ， 这 些 点 都 应 满足 恒等式 (22), 可 以 分 别 计算 出 每 个 点 的 值 代入 (22) 式 所 
得 的 相对 误差 . 随 着 n 的 增加 ， 近 奇 点 区 域内 的 点 也 会 增加 ， 相 对 误差 大 的 点 也 会 增多 .在 表 
2 中 列 出 了 相对 误差 的 统计 结果 ; 


表 2 IMSL 与 精细 积分 的 稳定 性 对 比 


点 的 个 数 相对 误差 > 10-12 的 点 的 个 数 最 大 相对 误差 

nxn IMSL 精细 积分 IMSL 精细 积分 
40401 8104 0 1.2910277291 x 10-? 4.8233240330 x 10-13 
161604 32380 37 7.7852746471 x 10-9 1.4551915228 x 10-11 
1010025 201709 228 2.3911707318 x 10-7 1.1641532269 x 10-1° 


从 表 2 不 难看 出 ， 精 细 积 分 法 无 论 是 精度 还 是 稳定 性 都 明显 的 优 于 IMSL 模块 . 
4 结论 与 展望 


本 文 利 用 精细 积分 法 对 二 阶 椭圆 函数 进行 了 成 功 的 求解 ， 并 通过 大 量 的 算 例 说 明 精 细 积 分 
法 具有 很 高 的 精度 与 稳定 性 . 但 目前 还 有 一 些 工作 需要 进一步 完成 : 

(1) Weierstrass 椭圆 滑 数 还 需要 作 更 多 的 误差 分 析 和 算法 对 比 ; 

(2) 将 精细 积分 法 应 用 于 其 它 特殊 函数 的 求解 . 
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